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Abstract 

We present results for the mass spectrum of cc mesons simulated on 
anisotropic lattices where the temporal spacing at is only half of the spa- 
tial spacing a s . The lattice QCD action is the Wilson gauge action plus the 
clover-improved Wilson fermion action. The two clover coefficients on an 
anisotropic lattice are estimated using mean links in Landau gauge. The bare 
velocity of light vt has been tuned to keep the anisotropic, heavy-quark Wil- 
son action relativistic. Local meson operators and three box sources are used 
in obtaining clear statistics for the lowest lying and first excited charmonium 
states of 1 5o, 3 Si, Pi, 3 Po and 3 Pi. The continuum limit is discussed by 
extrapolating from quenched simulations at four lattice spacings in the range 
0.1 - 0.3 fm. Results are compared with the observed values in nature and 
other lattice approaches. Finite volume effects and dispersion relations are 
checked. 
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I. INTRODUCTION 



Lattice Quantum Chromodynamics (QCD) opens a gateway to the study of non- 
perturbative phenomena in the strong interaction world. To explain (and in some cases 
predict) why the elementary particles are as heavy as they are is not only exciting, it is 
also unavoidable in the validation of QCD as the standard model of the strong interactions. 
Unfortunately the lattice simulations are by no means cheap, it is typical for a project to 
take months or years to finish on the present fastest supercomputers. 

The study of heavy quarks demands even more computing resources than that of light 
quarks, while heavy quarks may be more interesting. As standard lattice actions break 
down when the lattice spacing a > where niQ is the bare quark mass, the imposition of a 
fine enough lattice spacing makes the studies of heavy quarks completely out of reach given 
current computing power. To treat heavy quarks, special lattice actions must be designed. 
The two dominant approaches are non-relativistic lattice QCD (NRQCD) and the heavy 
relativistic or Fermilab approach, and there is the newer anisotropic relativistic approach 
used in our work. 

The NRQCD approach |T]fJ] attempts to describe an effective field theory at low energy. 
Essentially the action is expanded in powers of the lattice spacing a, as standard lattice 
actions are, and in powers of the heavy quark velocity v 2 . Practically the NRQCD method 
works well for the spin-independent bb system made of bottom quarks. However, continuum 
extrapolation is impossible in NRQCD as the non-relativistic expansion requires arriQ > 1. 
Also, to study spin splittings, higher order terms have to be added to the action. For 
the cc (charmonium) system, present evidence |3]|4j] suggests that the NRQCD approach 
breaks down. The spin splittings in the charmonium spectrum do not converge when higher 
relativistic corrections are added to the non-relativistic action, or when quantum corrections 
are switched from one prescription to another. 

The heavy-relativistic or Fermilab approach [|[] incorporates interactions from both the 
small- and large-mass limits. For heavy quarks, the lattice action can be interpreted in 
a non-relativistic light. Yet as m§a — > 0, the action conforms exactly to the standard 
Wilson action for light quarks. This is accomplished without any constraint on the value 
of am , in contrast to arriQ > 1 in NRQCD and arriQ -C 1 in Wilson action. The Fermilab 
approach connects both ends smoothly. Concretely, its lattice action up to 0(a 2 ) lattice 
errors is simply the standard clover-term improved Wilson action without imposing space- 
time exchange symmetry. The coefficients in front of the covariant derivatives and the clover 
term improvement now appear in two copies, a temporal one and a spatial one. The difficulty 
is that, to achieve the elimination of 0(a) lattice artifact for heavy quarks, these coefficients 
must all be all mass dependent. 

The anisotropic relativistic approach |||7j goes one step further. As in the heavy rela- 
tivistic approach, the space-time exchange symmetry is not imposed on the lattice action. 
The key difference is, here in the anisotropic relativistic approach, the lattice itself is dis- 
cretized differently along the temporal and spatial directions with the temporal spacing a t 
chosen finer than the spatial spacing a s . Since no further symmetry besides the space-time 
exchange symmetry is broken, the anisotropic action has the same terms as the heavy rel- 
ativistic action does. Defining the true anisotropy £ = ^ as the ratio of the spatial to the 
temporal spacing, we may consider the heavy relativistic approach as the £ = 1 special case 
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of the anisotropic relativistic approach. In both approaches, the relativity of the lattice 
action is restored by a required tuning of the bare parameters in the action. 

There is one obvious benefit of moving to an anisotropic lattice: heavy meson propa- 
gators often die out very fast and thus leave too few time slices which are useful for mass 
fitting. With a finer temporal lattice spacing, this problem may be cured at relatively low 
cost. Another equally important benefit is that, with e^mo < 1 on an anisotropic lat- 
tice, the mass dependence of the improvement coefficients can be expected to be weaker, 
or absent all together as is the case for some coefficients classically. Since the numerically 
determined clover coefficients deviate considerably from their perturbative estimates, their 
possible weak mass dependence may allow us to avoid a difficult, non-perturbative numerical 
determinations. 

To be self-contained, now we review the theoretical framework laid out in |||7]]. 

II. THE ANISOTROPIC WILSON QCD ACTION 

The anisotropic QCD action is the sum of the gauge action S G and the fermion action 

s ( = 4 + 4 • (i) 



A. The anisotropic gauge action Sq 



On an anisotropic lattice the gauge action becomes 



- Yl ReTr [l-P as ,(x)]+f £ReTr [1-P st (x)] 

SO x,s>sf x,s 



(2) 



where £ D is the bare anisotropy, which equals to the true anisotropy £ = ^ only at the 
classical level. Note that the anisotropic (3 is the geometric mean of the /3's along the 
temporal and spatial directions, thus it corresponds to a coarser spatial spacing and a finer 
temporal spacing than given by its isotropic equivalent of same value. Here a s and a t refer to 
the actual "physical" lattice spacings as computed by examining the propagation of physical 
particles over distances of many lattice units. Standard renormalization argument guarantee 
us that the resulting long distance physics predicted by the action in eq.(fj) will appear 
consistent with relativity after this anisotropic interpretation of lattice scales is adopted. 

It is suggested that the true anisotropy £ be fixed during the continuum extrapolation, 
for reasons shown later. All the computations are done at £ = 2 in this work. The bare 
anisotropy £ is tuned at each (3 to keep £ the same M. 



B. A few standard definitions 



The covariant first- and second-order lattice derivatives V M and A M are defined through 
their operations on the quark field q(x) 
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2a, 



U^(x)q(x + /J,)- U^ tl {x)q(x - /i) 



A>i*q( x ) = — U^(x)q(x + + U„^(x)q(x - jj) - 2q(x) 



Here we employ the notation U-^x) = U^x — jiy for the parallel transporter from x to 
x — 11. The lattice spacing 4- vector a M = {a tl a s l} is introduced to simplify the formulae. 
The Euclidean gamma matrices and the Dirac matrices a^ v are defined by 



T fiu — 2 [in ' lA ■ 



1p = %, {l^,lu} = 25^ , 
The field tensor F fiu (x) is defined by 

±Qnv(x) = U (x, n)U(x + jji, u)U\x + z>, n)U\x, v) + 

U (x, v)U\x — fi + z>, n)U\x — jl, v)U(x — jl, /J,) + 

U\x — jl, fx)U^(x — jl — v, i>)U(x — ji — v, fi)U(x — v, v) + 

W (x — v, u)U(x — v, fi)U(x + jl — v, v)U\x, fi) 



(3) 



(4) 



C. The anisotropic fermion action S F 



Back on an isotropic lattice, the following terms make up the clover improved quark 
action 



C 

1 il<V 



(5) 



in which the lattice spacing a is introduced to keep these terms of same dimension. 

For either an anisotropic lattice, or heavy quarks on an isotropic lattice, the space-time 
exchange symmetry should not be imposed at the level of the lattice action. Instead, in order 
to achieve a relativistic dispersion relation between energy and momentum, the coefficients 
in front of spatial terms in the fermion action have to be different from those in front of 
temporal terms. Thus, the anisotropic quark action is expected to be simply the standard 
isotropic action in duplicates, one spatial copy and one temporal copy, which is indeed 
exactly the final form we choose: 



Sf = a ta 3 s J2 l( x ) 



z s s<sr 



q(x) 



= ata 3 s Y <l( x ) 



m + z/ t A Wils ° n + E^P 



Wilson 

s 



S<s! 



q(x) . 



(6) 
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In the second line above, the first- and second-order lattice derivatives are combined into 
the r = 1 Wilson operators D™ lson = V M - ^f7 M A M . 

There are more bare parameters than in a standard quark action. The two clover coef- 
ficients are labelled by C* w and C s s w . The bare velocity of light v s is to be tuned to restore 
relativity on an anisotropic lattice, which could be equally well achieved by adjusting v t . 
Indeed we need to vary only one of them as the two cases are simply related by rescaling the 
quark fields. For later convenience we keep both v s and v t here. In practice we will always 
choose one of them to be 1 and tune the other quantity via the dispersion relation between 
meson energy and momentum. We will refer to these cases as "ivtuning" and "z/ t -tuning" . 

All the quantities in eq.(||) are dimensionful with hidden a s or a t . In order to program the 
quark action, we find it convenient to rewrite it in dimensionless quantities, i.e. quantities 

with a hat. With m = a t m , q = alq, V M = a M V M , A M = a jA„, D™ ™ = a^D™ *™ and 

F^y = ayfiivFfu,, the quark action reads 



4 = £g 

X 



x 



SO s A s 



s<st 



q(x). (7) 



Note that we choose to use £o instead of £ = — m the action. While £ = £ holds only 



classically, the another choice would just redefine v s and 



III. CLASSICAL IMPROVEMENT OF THE ANISOTROPIC ACTION 

The gauge action is already correct up to 0(a 2 ), so we only need to improve the quark 
action by choosing the right values for the bare parameters. Below are all the possible terms 
up to 0(a) in the quark matrix 

mo, yj t (l + O(am )), £ s ^(1 + O(am )), (8) 
aA t , aE s A s , oE^Jt], oE^Jt}, aY, s<sl {Ps,Ws,} ■ 

Note the term aJ2 s \fls,Wt] never arises on an isotropic light-quark action due to space-time 
exchange symmetry, and the last two anti-commutation terms are simply the clover terms 
in different faces. Also note the lattice spacing a has been reintroduced here to keep track 
of dimensions, referring to either a t or a s . All the coefficients in front of these terms are 
possibly mass-dependent. 



A. Field redefinition 

On the classical level the simplest way to derive the on-shell 0(a) improved anisotropic 
quark action is to relate it by a field redefinition to an action that has manifestly no 0(a) 
discretization errors. Since a field redefinition is just a change of variable in the path integral, 
on-shell quantities are not affected. The Jacobian of a field transformation matters only at 
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the quantum level, where, in the case at hand, its leading effect is to renormalize the gauge 
coupling and perhaps the bare anisotropy. 

We start with the naive fermion action which has no 0(a) discretization errors and whose 
bare quark mass is the continuum one m c 

q c (x)[m c +J7}q c (x) . (9) 

Then we apply the field redefinition q c = qQ, q c = qfl with 

fi = 1 + ^a t m c + ya t JJ t + ^-a$ (10) 

fi = i _|_ ^-a t m c + y a t JJ t + ya s ^ , 

where £l m> t, s an d & m ,t,s are six pure numbers, possibly mass- dependent. 

With all terms in eq.(^) showed up, the quark matrix in terms of the new fields q(x) and 
q(x) reads 

fi [m c + f] n (ii) 

/-, . fim fim \ , /i i i ^m fim . ^t\\ rn 

= (1 + a t m c ) m c + (1 + (a t m c h a t m c )) J/ t + 

/-, . / fim fim , fis fis \ \ m 

(1 + {a t m c h a s m c ))y + 

K^^)A, + (a^)(EA s + ER,W) + 

fi s — fi s fi t — ^ ^ fi s + fi s fit + fit ^ 
(a- 4 a < ^ ) ^ > ^] + ( a * 1 + a t ~ A ) \W , Vt] ■ 

A few words on the notation here: to avoid excessively complicated, rigorous expressions, 
we have simply pulled fi s , fi s and a s outside the spatial summation J2 S to make it explicit 
that they are common factors to all the three spatial directions. 



B. The classical 0(a) estimates of bare parameters 

To arrive at the anisotropic quark action in eq.(||), we first get rid of the expensive and 
meritless term [X 7 , X 7 *] i n ec L-(FD by requiring fi m = fi m , fi t = fi t and Q s = Cl s , which sets 
no limitation on the remaining terms. Now if we recall that the clover terms 

o^F^ = , J? v } , (12) 

the quark matrix eq. (|TT|) becomes 

fi [m c + X 7 ] fi = (1 + a t m c Q m ) m c + (1 + (a t m c Q m + a t m c Q t j) JJ t + (13) 
(1 + (a t m c Q m + a s m c Q s )) ^ JJ s + (a t fi t )A t + 

s 

(a s tt s )(J2A s + J2 Vss>F ssl ) + (a s ^- + a t y) ^a st F st . 

s s<sf ^ * s 
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We are left with three free parameters Q m , f2 s and Qf In the most attractive scheme, 
which actually leads to our definition of the quark action Sp in eq.(||), fl m is adjusted so 
that one (and only one) of u s or v t equals 1, Q s and Q t are adjusted so that the first- 
and second-order lattice derivatives are combined into the Wilson operators with the full 
projection property £)Wiison _ _ 

In the case of z/ s -tuning, the three parameters are set as Q m = |, Qt = — \ and fl s = 
~o !Ti° tmC ) giving the quark matrix as 

n [m c + W n = (1 + l -a t m c ) m c + P t Wilson + i±|^£ £ pWiison _ (M) 
2 1 + 2a s m c V 

a ss/F ss/ - -(-^ + -^) • 

Z s<s/ Z Z Z s 

From eq.(|14D we can read off the classical estimates of the bare parameters in the anisotropic 
quark action eq.@ if that action is to have no 0(a) errors 

mo = m c (l + \a t m c l v s = \ + »t = 1, C* w = 1, = hi + ^) . (15) 

£j ±- I 2 s c s 

In the case of ^-tuning, the three parameters are set as Q m = f2 s = — | and 
fit = — ~ ^ + 2 as ' Ttc _ ^he c l aS si ca l estimates of the bare parameters are 

m = m c (l + \a s m c ), v s = 1, v t = \ + f' m \ Q w = 1, <£, = \(l + ^) . (16) 
2 1 + ^a t m c 2 a s 

The parameters specified by eq.fllSD and eq.fllTf) correspond to our two possible con- 
ventions for v s and v t . Since these two conventions are connected through a simple scale 



transformation of the fermion field q(x) — > \J Vt/v s q(x), the parameters in eq. (|I6"D should 
equal those in eq.(fl5|) multiplied by ^ = ^f~~ ■ Of course, this simple scaling holds only 
to 0(a), a limitation that will be improved below for the clover coefficients C* w and 



C. Better estimates of the clover coefficients 



Since the 0(a) dependence of the clover coefficients leads to only 0(a 2 ) errors of the 
action and we only aim to remove the 0(a) lattice artifact, when writing down the classical 
estimates of C* w and C s s w in eq . (|14|)- (|16|) we have neglected the 0(a) parts of the transfor- 
mation coefficients Q s and Q t - However, the neglect of 0(a 2 ) terms in this manner is not 
necessary, indeed the clover coefficients are better expressed in terms of the bare velocity of 
light v t or v s . In this way, we can partially determine the mass dependence of the clover 
coefficients and also resolve the contradiction that the clover coefficients given in eq . (|15"D 
and eq.flTB]) are not related by the factor ^. 

This statement is based on the following observation on the general form of the quark 
action in eq.([L3|). We see that the spatial clover term Y, s < s / a ssiFssi always comes together 



7 



with the spatial Wilson term J2 S A s and there is also a similar relation between the temporal 
ones. Therefore, with the same values of Q m , Q m and Qt given in above paragraphs, we can 
give a more precise version of eq. (p!5|)-(p^). In the ^-tuning case, the more precise classical 
estimates of the bare parameters 

m = m c (l + \a t m c ), v s = a * 771 " , v% = i ; C* w = u s , C* w = ~(^ s + — ) . (17) 

In the ^-tuning case, the estimates are 

m = m c (l + -a s m c ), v s = 1, v t = \ s \ C s sw = 1, C* w = -(1 + u t —) . (18) 
l l + -^a,tfn c i a s 



IV. COMPUTATIONAL PROCEDURE 

We have run quenched simulations at four values of the lattice spacings (see the input 
parameters listed in table [I]). The mass spectrum in the continuum limit is then obtained by 
extrapolating the measurements at these finite lattice spacings to zero lattice spacing (see 
results in table [V]). 

A. Adapt existing isotropic software to anisotropy 

Fortunately it is trivial to modify existing isotropic software to simulate anisotropic 
lattices, at least for quenched calculations. Our practice is to re-scale the temporal links so 
that the anisotropic lattice action appears like the standard isotropic action. For the gauge 
sector, the temporal links are multiplied by £o> transforming the gauge action from eq.(^) 
into 

4 Ut ^= Ut tw £ ReTr t 1 - p ^( x )\ + const • ( 19 ) 

For the fermion sector, the temporal links are multiplied by ^r£o> transforming the quark 
action from eq. (|7|) into 



m + f- p Wilson - (20) 



1 v 2 C s 

/ u t<,0 s ?0 s <sr 



q(x). 



In this way, we can use existing heat-bath code to update links and existing Dirac operator 
^9 Wilson to measure mass spectrum, as long as the code does not assume the SU (3) properties 
of the gauge links, for example, to reconstruct third row of these links. The two different 
scalings are not a problem in quenched simulations, although it may require some thought 
in full simulations as both Sq and Sp are used in combination to update field configurations. 
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It is a popular practice in lattice simulations to write the quark action in terms of k and 
thus separate local terms A from nonlocal terms p^ lson 



PT° n q(x) = 7T- Et(7. - r)U,(x)q(x + - ( 7 „ + r)U.,(x)q(x - /x)] 



(21) 



On anisotropic lattices, a/ter the rescaling of temporal links U t — > (,o^U t , the quark action 
can be rewritten as 



Si 



u t 1 



. Wilson, 



(22) 



where the anisotropic k and A are defined as 

1 



K 



A* 



2 

1 - i; 



, in this work r = 1 



7/, 



(23) 
(24) 



B. Set the lattice scale 

The heat bath algorithm for updating the gauge configurations is the standard Creutz 
method extended by Kennedy- Pendleton [10] and Cabibbo-Marinari [llj]. Moving from 



an isotropic lattice to an anisotropic lattice, we need two bare parameters, namely the bare 
anisotropy £ an d P, to specify the gauge couplings needed to generate the gauge links. In 
order to extrapolate the measurements to the continuum limit and to run simulations at 
reasonable physical volume, we need to know the lattice spacings a s and a t in physical units, 
or alternatively the spatial spacing a s and the renormalized (true) anisotropy £ = a s /a t . 
This work has been done in || and fll2|-(14}| . The part of these earlier results used in our 
simulation are listed in table 0. We now briefly describe their methods. 

In our simulation the renormalized (true) anisotropy £ has been fixed at £ = 2. The 
choice of £ to be an integer makes it most convenient to determine the relationship between 
£ and £ at a given P ||. Basically two static quark potentials are compared with each 
other. Both quark-antiquark pairs are propagating in the same spatial direction. One pair 
is separated in a spatial direction (different from the propagating direction, of course) while 
the other is separated by twice (or other value of £) the number of lattice sites in the 
temporal direction. The bare anisotropy £o is tuned until the two static quark potentials 
become identical. 

It is also a good idea to be cautious by keeping the value of the renormalized anisotropy 
£ fixed during the course of taking the continuum limit. In this way the scaling property 
of the mass spectrum is the only issue here and we avoid the complication all together that 
lattice physics might behave differently at different anisotropy. Beyond this reason one can 
always take the continuum limit along some other smooth curve of true anisotropy £ varying 
with lattice spacing a. 
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The lattice spacing was determined very accurately in terms of the Sommer scale r$ 
14 . The scale r$ is defined via the force between a heavy quark and antiquark 



r 2 Q F(r ) = 1.65 . (25) 

The constant 1.65 is chosen so that r = 0.50 fm is an intermediate distance quantity. By con- 
trast another often used dimensionful quantity, the string tension a, is only asymptotically 
defined at long distance and thus suffers from inexact assumption of leading intermediate 
distance corrections. Although both the Sommer scale and the string tension a are cal- 
culated from the static quark potential, the available data are much better at intermediate 
distance. Therefore it is superior to use = 0.50 fm to attach physical units to lattice 
observables. 

We should note that although all successful potential models closely agree on the value 
of ro to be 0.50 fm, it does bring some systematic errors to our simulation. We decide not to 
quote any unjustified estimate of the errors on r = 0.50 fm and leave it open to the reader 
as to how to treat the errors from this choice of length scale and the effect of quenching. 
We should also note that this error comes in only at the very end of our calculations when 
the mass spectrum extrapolated to zero lattice spacing is written in the physical unit of 
MeV. Note, both ^ and £ = £o(£>/3) were determined to 1% accuracy. Therefore when 
extrapolating in the lattice spacing we may neglect the errors from ^ and £o and worry 
about the errors on our mass measurements only. 



C. Estimate the clover term coefficients 



On each gauge configuration generated with the heat-bath algorithm, we invert the 



fermionic matrix using the conjugate gradient (CG) method [15| with preconditioning. We 
then measure all mesonic states that can be obtained from bilinear sources without derivative 
operators (using the later would lead to more noisy correlators), as shown in table |T|. 

In addition to the bare quark mass, these measurements require three more input param- 
eters in the fermionic action. Two of them are the temporal and spatial clover coefficients 
C* w and C s s w . Since a non-perturbative determination of the clover coefficients (say using 
the Schrodinger functional) is a daunting project, C* w and C s s w are estimated using tree- 
level tadpole improvement. There is empirical evidence that tree-level tadpole improvement 
achieves more than two-loop or even three-loop perturbative improvement does. At tree- 
level tadpole improvement one starts with the classical action and then renormalizes each 
gauge link by its "mean value" 

U, - ^ . (26) 

Clearly the mean links = {ut, u s ,u s ,u s } can not be defined in a gauge invariant manner. 
The prescription to isolate the true gauge-independent tadpole contribution is to minimally 
renormalize the gauge links by choosing a maximum definition of the mean links u^. 

This maximization of mean links leads us to determine them in Landau gauge since on 
the lattice, the Landau gauge condition d^A^ = is achieved by maximizing the functional 
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a 2 



(27) 



However, there is one subtlety regarding the ratio a t /a s of spatial and temporal lattice spac- 
ing: which anisotropy should be used in this gauge fixing process, the bare or renormalized 
one? We choose the bare one £o based on the following empirical observation P|||7|]. The 
tadpole improvement hypothesis says that the ratio Ut/u s gives a tree- level estimate of the 
renormalization of the anisotropy £ = — £o (this relation will be used in the next paragraph 
to simplify eq. (|29|) ). It is found in |||7j that with the choice of £o the measured u t /u s in Lan- 
dau gauge agrees quite well (within 2% or less) with the values measured non-perturbatively 
in!- 

The tree-level estimate of the clover coefficients is given in eq.flTSD-flnSD as 



a 



t 

sw 



2 V + l 



(28) 



Now we work on the tadpole correction starting from the quark action 



tadpole 



m + v t A Wilson 



S s 



T^swE "*^ + E a ssfF ss/ ] 



s<sr 



q(x). 



Et( 



X 



rn Q + L Vt ^Wilson + ly_L y ^Wilson. 
Ut U s £ V 



2 uiu s 



u 



1 n s 

/ , Vssi r ssi\ 



s<s/ 



q{x) 



Eft 



w t m + i/ t p t 



Wilson 



6 



E£ 



Wilson 



11 1 c s 

q [ J^sw E + ~~3 "T^ E °"ss/-Fss/] 

2 u t wf V uj £ 



s<s/ 



g(x) 



(29) 



Comparing above with the action eq.(0) we simulate, the tree-level tadpole estimate of the 
clover coefficients is 



l 



a 



w 



t 

sw 



2 i 1 + J ttt«2 



(30) 



The values given in eq.(|30D are what we use in this work. However, in retrospect, we think 
it may be better not to base the tadpole improvement on eq.(|ll^-(|T6|) or eq. (|28|) . Rather, 
there is a better classical estimate given in eq.([l7)-(|l8). In the ^-tuning, the difference 
between eq. (|15D- (pT6|) and eq.([I7|)-([l8[) is small and perhaps negligible, but this doesn't seem 
to be the case for the ^-tuning. We will comment more on this in section |V H. 
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D. Tune the quark mass m,Q and the bare velocity of light u t 



The remaining bare parameters in the quark action are the bare velocity of light v t and 
the bare quark mass m . Before we explain how these two inputs are tuned, we need to define 
the effective velocity of light c(p) first. In terms of the energy E(p), the mass m = E(0) 
and the momentum p of a meson, the dispersion relation has the form of 

E\p) = m 2 + p 2 + 0(£vt) + 0(ZpIpD- (31) 

s ss' 

The effective velocity of light c(p) is then given by 

c(p) = t 



a 2 E 2 (p) - a 2 m 2 



where the factor £ = a s /a t comes from the fact that the lattice energy is expressed in 
temporal spacing a t while the momentum p is expressed in spatial spacing a s . The fact that 
c(p) 7^ 1 is due to finite lattice spacing error. 

The bare quark mass and the bare velocity of light u t are tuned simultaneously so 
that the spin average IS meson mass equals its observed value 

im(l 3 S ) + |j™(l 3 Si) = 3.067 GeV (33) 

and c(0) = 1 for the pseudo-scalar meson rj c . 

We obtain c(0) by extrapolation from the c(p) of the pseudo-scalar meson r\ c at the two 
lowest, on-axis momenta ^(1,0,0) and ^(2,0,0), assuming c(p) — c(0) oc p 2 (we may 
also choose some other direction in momentum space, say ^(1,1, 0) and ^-(2, 2, 0), but the 
0(p 4 ) errors in c(p) will be larger). To increase statistics we always average over momenta 
that just differ by permutations of their components. Precisely we use 



c(0) 



16£(1) 2 - £(2) 2 - 15m 2 
\ 12(2tt/L) 2 ' (34) 



where the three energies E(2) = E(p = ^(2,0,0)), E(l) = E(p = ^(1,0,0)), and m = 
E{p = ^-(0,0,0)) are from a correlated fit of hadronic correlators at three momenta. But 
this formula is merely one way of computing c(p). The point is we eliminate the leading 
relativistic errors completely by demanding c(0) = 1 up to 0(p 4 ) error. 

The tuning here involves two to four iterations to get both the quark mass tuq and the 
bare velocity of light u t right to about 1% (see table |), an accuracy in line with that of 
scale setting a s /r 0N and £ - Generally, the simultaneous tuning of more than one parameter 
to such a precision can be quite expensive, but here is relatively easy, since we find that 
the mass dependence of v t is very weak on anisotropic lattices, just as in the classical case. 
Besides it is sufficient to use about 100 configurations for the initial tuning. For the final 
measurement runs using the tuned parameters we have over 400 configurations per lattice 
spacing. 
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In our experience, an increase of the bare velocity of light v t leads to smaller meson 
masses and effective velocity of light c(0); not surprisingly, an increase of the bare quark 
mass m leads to larger meson masses although such an increase has only very minor effect 
on c(0). Thus, one may tune the bare velocity of light v t first until c(0) = 1 and then work 
on the bare quark mass mo to make the meson mass m(lS) right. 

Alternatively we may set v t = 1 and tune u s instead. The best value of u s should be 
around the inverse of the best value of u t . The desired C* w , C^ w and m of z/ s -tuning are 
simply their corresponding value of f t -tuning divided by u t since two cases are related by 
rescaling the fields by a factor of u t . 

All together in this work we have studied the dispersion relation using the six lowest 
momenta. Omitting the common factor 2zc, they are 

pee = [0, 0, 0], [0, 0, 1], [0, 0, 2], [0, 1, 1], [0, 2, 2], [1, 1, 1], [2, 2, 2]. (35) 

While P3 and p4 are useful in double checking the tuning from P0..2, Ps and p 6 are too noisy 
to be of any use. 

E. Fit hadronic correlators of local sink and three box sources 

Our computation gives a very extensive mass spectrum, namely, the radial n = 1 ground 
states and n = 2 first excitations of the 1 Sq, 3 Si, l P%, 3 Po and 3 Pi particles as listed in table 
An anisotropic lattice certainly gives the benefit of fine temporal spacing without much 
computing expense. Had the simulation been performed on an isotropic lattice of equivalent 
cost, the signal of the heavy hadronic correlators may have died out quickly, making the 
mass fitting technically impossible. 

On each gauge configuration we compute the quark propagator three times. Each time 
only the size of the box source differs, as detailed in table [TV]. Each calculation uses local 
sinks. Although this combination is largely due to software availability, it is sufficient to 
generate masses of the ground state and first excited state with better precision than given 
by typical lattice calculations (see results in fig [l] and table [V|) . The trick is to gear the size 
of the box source to correspond to the desired wave-function size. Fortunately we need to 
do this tuning at only one value of lattice spacing. Since we know the lattice spacings in 
terms of physical units, we can then easily estimate the optimal box sizes for other values 
of lattice spacings. 

In theory a fitting ansatz should incorporate the energies -E'i(p), E 2 (p), etc of each state 
entering the hadronic correlator with point source y and point sink x 

< ]T e^ p - x ft(x, r)H(y, 0) > = £ | < n\H\0 > \ 2 e~ En{p)T . (36) 

x n 

However, we only have a limited number of time slices and therefore we want to reduce the 
number of fitting parameters as much as possible, so long as the fitting ansatz still closely 
reflects the time dependence of the underlying hadronic correlator. 

The three sizes of box sources at each lattice spacing will be referred as small-size, 
medium-size and large-size. A 1-cosh fitting ansatz (i.e. ground state only) applies to the 
hadronic correlators with the medium-size box source so well that the fitted ground state 
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mass stabilizes for a fitting range as early as the minimum time slice t m i n = |r. The other 
two hadronic correlators are to be fitted with a 2-cosh ansatz or a 3-cosh ansatz to give 
masses of excited states. We did not use hadron correlators with a point source because 
the undesirable contributions they receive from higher excited states (of radial quantum 
number n > 3) can not be overlooked. Instead, we speculate that a small-size box source 
behaves like a "mild" point source. Just as | < n|7Y|0 > | 2 in the case of a point source, the 
contributions of each energy state in the hadronic correlator of a small box source are all 
positive. Nevertheless, the correlators containing a small-size box source do not appear to 
be contaminated by higher excited states beyond our interest. Meanwhile the correlators of 
a large-size box source may behave like a wall source correlator. The amplitudes of different 
energy states may be a combination of positive and negative numbers. When fitting three 
correlators simultaneously it is indeed a desirable feature that a physical state manifests 
itself in these correlators with different signs. 

Powell's method is used to minimize the x 2 to generate the best-fit parameters. Devel- 
oped by Kent Hornbostel and Peter Lepage, the fitting code supports an arbitrary (corre- 
lated) fitting ansatz on arbitrary numbers of data files. We use 1-cosh, 2-cosh and 3-cosh 
fitting ansatze. The t max of a fitting range is typically fixed through the effective mass cal- 
culation while the t m \ n is varied through all values as long as there are enough degrees of 
freedom left. 

To be selected, a fitted result ought to satisfy three criteria: 1. It is consistent with 
all other fitting ansatze. For example, although the 2-cosh ansatz reaches plateau earlier 
than the 1-cosh ansatz does and only one of them may succeed in generating meaningful 
best-fit parameters, the ground state mass should agree between both ansatz. 2. It becomes 
stable when t min reaches certain value. 3. All the fitted quantities including amplitudes are 
statistically nonzero and have the right sign if known. Precisely the fitting procedure is: 

1. Calculate the effective mass on each adjacent pair of time slices using the 1-cosh ansatz. 
This step supplies the value of t max and the estimates of ground state masses and am- 
plitudes to step 2. It also gives us an opportunity to easily examine the autocorrelation 
between configurations. 

2. Apply the 1-cosh fitting ansatz on each individual datafile to give better initial guesses 
for the real fittings to follow. One datafile corresponds to one unique combination of 
mesonic operator ipTip, momentum p, and box source size. 

3. The spin average mass m(lS) and the mass differences Am^s^i^, Am^^, and 
Am 1 3 Pl _ 1 3p are obtained using the correlated 1-cosh ansatz on two or three (i.e. the 
number of particles involved) correlators at zero momentum, p = 0, and medium-size 
box sources. Only the correlators of medium-size box sources are used here since they 
are designed to suit the 1-cosh ansatz very well. If we otherwise include correlators of 
small-size box sources and large-size box sources, we will have to add more states into 
the fitting ansatz. The fitting will then become too complex to succeed. Since one of 
the triplet P-wave states, 1 3 P2, is missing, we do not have the true spin-average of the 
l 3 Pj states. Instead we have used l l P\ in the l l P\ — IS splitting. 

4. To monitor the effective velocity of light c(p = 0), a correlated 1-cosh ansatz is ap- 
plied on three 1 S , correlators of momentum 0, ^(1,0,0) and ^(2,0,0). Again, only 
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correlators of medium-size box sources are used here. We observed however, that the 
higher the underlying momentum is, the smaller would be the best medium-size source 
for a 1-cosh ansatz. Whatever size we choose for the box source, it will not make the 
correlators of all momenta simultaneously perfect for the 1-cosh fitting. But since the 
effect is weak and the fitted results are stable across different t m i n , we should not worry 
too much about the relatively worse x 2 /d.o.f. here. 

5. To obtain the masses of the ground state and first excitation of each particle 2S+1 Lj, 
the correlated 2-cosh ansatz is applied on its correlators from small- and large-size box 
source while a 1-cosh ansatz is used simultaneously for the medium-size box source. We 
always apply 1-cosh ansatz on the correlators with medium-size box source because the 
excited state amplitude of correlators from medium-size box sources are statistically 
zero when fitted with 2-cosh ansatz. 

6. The 3-cosh fitting is done for the purpose of a sanity check. 

The masses and mass differences are quoted in table [V] and [VI|. The fitting details 
such as xVd.o.f. are listed in table |VH| - |XI| . In these tables the Q value is a normalized 



indicator of the quality of a fit, defined as the probability that we would end up with a 
higher x 2 (Xmin ^ more strictly speaking) if we did the simulations many times, it takes 
values between and 1. The dropped eigenvalues refer to the truncated smallest eigenvalues 
of the sample correlation matrix. The goodness-of-fit is chosen as the product of the Q value 
and the degrees of freedom (d.o.f.). The d.o.f. in turn is the number of time slices from all 
correlators minus the number of fitting parameters and the number of dropped eigenvalues 
of the correlation matrix. 

The error of a fitted quantity is given as the amount of perturbation away from the 
best-fitted value in order to increase the \ 2 by 1. Assuming the data model is right, this 
definition indeed gives the 68% range of a fitted parameter [IBj. And it is much faster than 



the bootstrap or jackknife methods because it avoids the labor of producing and fitting 
synthetic data sets. For a poor fit, usually signaled by a large x 2 /d.o.f. > 1.5, this definition 
may underestimate the statistical error. The majority of our fits have x 2 /d.o.f. ~ 1, thereby 
we expect the errors from A% 2 = 1 are good estimates. It is also a good idea to look at the 
fluctuations of the fitted results under a variety of fitting conditions and then to take these 
fluctuations into account in quoting the statistical errors. 



F. Extrapolate to the continuum limit 

Ultimately we want to reach the continuum limit by extrapolating from the computations 
at four values of lattice spacings. Only then can we compare our results with, and in some 
cases predict, the mass spectrum in Nature. 

The mass differences, not the masses themselves, are to be extrapolated to the continuum 
limit for several reasons: 1. We have tuned the spin average meson mass m(lS) to equal 
their experimental value, therefore the masses will more or less right by design. 2. The 
tuning is not perfect, but instead is accurate to 1%. Consequently the mass spectrum from 
the computation at one lattice spacing may all be systematically lifted up, while from a 
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different run they may all be dragged down. Such overall mass shifts cancel out in the mass 
differences. 3. The masses are not independent quantities since they are all measured on 
the same gauge background. A correlated fit for the mass difference may be more precise 
than a naive subtraction of two masses fitted independently. 

The computations at different lattice spacings are independent computations, which 
makes the continuum extrapolation technically very easy. Both the Wilson gauge action 
and the clover-term improved Wilson fermion action are accurate to 0(a 2 ). Therefore the 
measurements at four values of a 2 s are fitted with a straight line (see table [V| and figure 
|2|-^p. At (3 — 5.7, the run on a 16 3 ■ 64 lattice is used due to its better fitting quality than 
that of the run on a 8 3 • 32 lattice. The intercept at a 2 s = gives the continuum extrapolated 
value. The extrapolation is done using the regression functionality of the software Xmgr, 
and is double checked using Maple. 

V. RESULTS FOR THE CHARMONIUM SPECTRUM 

Now we present our results and compare them with (if available) the experimental data 
and other lattice approaches. 

A. Finite volume effects 

In this work we have run simulations on four values of lattice spacings. As is important 
in every lattice calculation, we need to make sure the simulation volumes are large enough 
to avoid significant finite-volume effects yet also small enough to avoid high computational 
cost. Among the four (3 values, the (3 = 6.1 run has the finest lattice spacing, its lattice of 
16 3 ■ 64 sites extends 1.536 fm along the spatial directions. Therefore we choose to check 
for the finite volume effects at (3 = 5.7 on two lattices 8 3 • 32 and 16 3 ■ 64, corresponding to 
a spatial extension of 1.658 fm and 3.315 fm respectively (see table 0). As listed in tables 
| and |I^, the identical input parameters are chosen for the two (3 = 5.7 runs. The mean 
links in Landau Gauge are measured on a lattice of 16 3 ■ 48 for the purpose of estimating 
the clover term coefficients. 

Comparing the mass spectrum between these two volumes at (3 = 5.7 (see tables |V| and 
|VlD , we see no sign of an overall volume bias and for the radial n = 1 ground states we see 
broad agreement within one sigma. The fitting of the 16 3 • 64 run is easier and more stable 
(over ttmin) than the fitting of the 8 3 • 32 run. The reason is unclear — presumably at the 
larger volume each hadron correlator fluctuates less after the time slice average. However, 
we do expect both fittings to work better if the smallest box source size 2-2-2 was slightly 
larger, say 2-3-3. 

The two c(0)'s are consistent within 1% (see table |). Note the bare velocity of light v t 
is only tuned at the smaller volume 8 3 • 32 by demanding c(0) = 1, this value of v t is then 
adopted in the simulation at the larger volume 16 3 • 64. Also note that c(0) is extrapolated 
from the energies at momenta 0, ^(1, 0, 0) and ^(2, 0, 0). Since the spatial extension L of 
the 16 3 -64 run is twice that of the 8 3 -32 run given the same lattice spacing, the extrapolation 
of c(0) at the larger volume actually uses a different set of momenta. Thus the consistency 
between the two c(0)'s also serves as an important check of the dispersion relation. 
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The 16 3 • 64 run is used in the continuum extrapolation because its mass splittings are 
more accurate. 



B. Dispersion relation 

Compared with other lattice approaches to heavy quark systems, our approach has two 
distinctions: 1. the temporal lattice spacing is finer than the spatial lattice spacing; 2. the 
relativity of the lattice action (broken in all heavy quark approaches) is restored numerically. 
Recall that we tune the bare velocity of light u t so that the c(0) of the pseudo-scalar l 1 ^ 
state is close to 1. For the validation of our method of restoring relativity on an anisotropic 
lattice, it is important to check the universality of c(0). We have checked it once in the study 
of finite volume effects by comparing two c(0)'s on 8 3 • 32 and 16 3 • 64 lattices at /3 — 5.7. 
A more extensive check on the dispersion relation was done in the [3 — 6.1 run on a 16 3 • 64 
lattice. 

The universality of c(0) indeed holds (see table |X11|) for a variety of particles l 1 ^, 



1 Si and 1 Pq and holds for the extrapolations from two sets of momenta, one set is 
[0,^(1, 0,0), ^(2, 0,0)], another set is [0, ^(1, 1, 0), ^(2, 2, 0)]. Other particles and mo- 
menta are missing purely due to failures in data fitting. 

Compared with the "non-relativistic reinterpretation" (to be described in this para- 
graph), our approach of maintaining the relativity of an anisotropic or heavy-quark Wilson 
action is conceptually clean. Furthermore, it has a statistical advantage as well. It is argued 
that without restoring relativity, the dispersion relation may be expanded in velocity as 

E(p) = m static + p2 + 0(p 4 ) (37) 

kinetic 
P 2 

— m static ~r 



2m kine t ic (p) 

in which the kinetic mass mkinetic is the true mass m, yet a difference between two static 
masses Am static also provides the true mass splitting as does Am kinetic . The mimetic is 
obtained in the same way we use to obtain c(0), namely by extrapolation from c(p) or 
^kinetic (p) at the two lowest on-axis momenta. Combining the c(p) equation 

£(p) 2 = m 2 + p 2 c( p ) 2 (38) 

with the eq.(|3~7D gives 

^■static /Qn\ 

m kinetic = C ^Q^2 • l""J 

As checked in ||[7j], eg . (|39|) is indeed true within errors and the relative errors of mimetic 



are indeed roughly twice those of c(0) as predicted by this equation. The absolute errors 
of 

^kinetic are a lso found []6],|7j to be one order of magnitude larger than those of wi s tatic- 
This difference in statistical errors is not hard to understand: the kinetic mass comes out 
of a correlated fit on three meson correlators (i.e. from three momentum values), while the 
static mass is given by the meson correlator for zero momentum alone. Presumably in the 
calculation of hadronic correlators, there is also extra difficulty in finding a smearing or box 
source size that is good for all three momenta. 
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C. The charmonium spectrum 



In figure [I] we plot the mass spectrum from quenched simulations on four anisotropic 
lattices against the experimental values [|17]], and in table |V| we list the precise numbers. We 
can see that at this scale the agreement between our lattice simulations and the observed 
values in Nature is very impressive and that even the effect of quenching is hard to see. 

The masses come out of the computer in units the inverse of the temporal lattice 
spacing. Since the lattice spacing is never an input and is always treated as 1 in simulations, 
in order to quote the masses in the physical unit GeV or MeV, we have to know the physical 
size of the lattice spacing at (or equivalently that of a s as the true anisotropy £ has been 
fixed at 2). In addition, the bare quark mass has to be tuned to correspond to the charm 
quark mass so to obtain the charmonium mass spectrum. In determining the lattice scale a 
we have not used any meson masses, instead, we set the lattice scale using the Sommer scale 
r because it can be measured more accurately than the popular choice of the l l P\ — IS 
mass splitting. To fix the bare quark mass, we set the spin average IS" meson mass m(lS) 
to its experimental value. All the remaining energies, both the n = 1 ground states and the 
n = 2 excited states, are then predictions. 

Among these predictions, most of the ground states deviate from experiment by less 
than 30 MeV. Regarding the minor discrepancies seen among the four estimates of each 
ground state, at least 50% may be attributed to the initial choice of quark mass (which is 
tuned accurately to 1%). For example, there is a downward shift on all the masses from the 
(3 = 6.1 run because the initial quark mass is slightly too small. For the excited states the 
deviations from experiment are typically under 100 MeV. Note there are no experimental 
data available for the excited states of particle h c , Xco and Xci- 

The statistical errors on the excited states are one order of magnitude larger than the 
errors on the ground states for the following reasons. The signal of an excited state, being 
proportional to exp(— m excit ed), dies out much faster than the signal of a ground state, there- 
fore far fewer time slices are useful in a fit determing the excited state mass. Furthermore, 
in our calculations the excited state signals are always mixed with the larger ground state 
signals, making the fitting of an excited state subject to the errors from the fitting of the 
ground state. For reasons mentioned earlier, the errors listed in the data tables are purely 
statistical errors, thus they do not include the systematic errors from the Sommer scale 
setting, the quark mass tuning, or from quenching. 



D. l 1 ^ - IS splitting 

In lattice simulations of heavy quarks the l*Pi — IS splitting is often used to set the 
lattice scale, denoted as a-iip^s- Because the Sommer scale r can be measured more 
accurately, We have used r = 0.5 fm to set the scale, denoted as a ro . In order to see how 



these two methods of scale setting differ, we plot our results (see figure £| and table |VI[) in 
the form of their ratio a 1 ip 1 _ 1 5/a ro 

a^Pi-is _ Amii Pl _ig 

a ro ~ 458.5MeV 1 ' 
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where Ampp^^ is the l x Pi — IS splitting with physical units from setting ro = 0.5 fm, 
and 458.5 MeV is the experimental value for the l x Pi — IS splitting. 

The continuum a 2 extrapolation gives the ratio an^-is /Or — 0.94(1). This discrepancy 
from 1 may come for two reasons. One is due to quenching, the splitting l l P\ — IS is 
smaller than its physical value, therefore aiip 1 _i5 has been underestimated. Another reason 
is associated with ro = 0.5 fm, which is only a phenomenological estimate and not a hard 
experimental number. Any errors with the assignment ro = 0.5 fm will affect our final values 
of masses quoted in physical units but only these final numbers. We also show the most 
recent results based on the Fermilab approach (|I£,nj for comparison, where the agreements 
are obvious. 



E. l 3 5i - l l So splitting 

Here comes the most exciting part of our results: the hyper-fine structure of the char- 
monium mass spectrum, plotted in figure |5| and listed in table [VI]. From the continuum a 2 
extrapolation, the mass splitting Am 1 3 Sl _ii So comes out to be 71.8(20) MeV, which is 39% 
smaller than its observed value of 117.1(2) MeV in nature. 

For really heavy quarks, one may speculate that the hyper-fine splitting will be dominated 
by one-gluon exchange. Thus one might be able to estimate the quenching effect by looking 
into the difference in the running of the strong coupling constant in quenched and full QCD. 
The 39% discrepancy from our simulations supports the expectation of large quenching 
effect on the hyper-fine splitting, although it is probably too aggressive for us to claim a 
39% quenching effect without qualification since we have not numerically tuned the clover 
term coefficients. 

As to the comparison with other lattice approaches, our results are consistent with 
calculations from the Fermilab approach fl8|]l9|| . Also shown in figure [5] are the NRQCD 
results at their best. As realized and fully discussed in HU, the NRQCD results can 
not be trusted. Inconsistent results are given by actions which differ only in the order of 
relativistic correction or which differ only in the tadpole prescription for quantum correction. 
Therefore the NRQCD results will not be included in our later comparisons. 



F. l 3 Pi - 1 3 P splitting 

The P-wave fine structure of the charmonium mass spectrum is shown in figure with 
exact numbers listed in table [VT|. The continuum extrapolation says the mass splitting of 
l 3 Pi — l 3 Po is 65(3) MeV, which is 30(5)% smaller than the experimental value of 93(3) MeV. 
As in the case of >S-wave splitting, most of the discrepancy with experiment is attributed to 
the quenching effect and it is not hard to claim consistency with results from the Fermilab 
approach [ p~8| , |1^ 1 . 
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G. Effects on results from small changes of bare parameters 



Here we discuss how the outputs (masses, mass splittings and the relativity indicator 
c(0)) respond to a 5% or 10% change of simulation inputs. We will examine four inputs: 
the bare quark mass m , the bare velocity of light v t , and the two clover term coefficients 
Cg W and C* w . In our calculations we have tuned the inputs v t and mo numerically so that 
c(0) = 1 and m(lS) = 3.0676 GeV. By contrast, inputs C s s w and C* w are estimated (not 
tuned) from mean links in Landau gauge using tree- level tadpole improvement. To help the 
tuning of v t and mo and to estimate the effects from the absence of numerically determined 
Cg W and C* w , we need to know the quantitative sensitivity of our results to these inputs. 

Detailed in table |XIII| and identified as run to 6, seven tests are done for this purpose. 
Run is to be compared with others. As to the remaining six, each of them differs from 
run only in one input parameter. Now let's look at the results in table |XIII| to see the 
effects of each input one by one. We have listed three types of mass splittings in the table, 
namely the spin-spin splitting Am 1 3 Sl _ 1 i So , the spin-orbital splitting Am^.^, and the 
S — P splitting Am 1 ip 1 _ 1 £'. However we will only focus on the spin-spin splitting as the 
characteristics of the other two are either the same or hard to tell given their relatively larger 
errors. 

In runs 1 and 2, the bare quark mass m is changed by ±5% from its best value of 0.51 
as tuned and used in run 0. The comparison of these three runs shows that mo has no effect 
on c(0) or on mass splittings. However, as expected, an increase of the bare quark mass mo 
gives a boost to masses of its bound states, m(lS) listed as an example. The errors from 
the tuning of mo therefore drop out of the discussions of mass splittings. 

In runs 3 and 4, the bare velocity of light v t is changed by ±5% from its tuned best 
value of 1.01. Comparing run 0, 3 and 4, we see that an increase of v t reduces the values 
of masses, mass splittings and c(0). This observation agrees with what is indicated through 
a field redefinition: putting the bare relativity factor in its more conventional place, i.e. in 
front of the spatial derivative, 



q(x) [mo + v t Pt+ P] q{x) — ► q{x) 



mo ^ 1 ' 

— + A + - P 

vt v t 



q{x) (41) 



we see that effectively ^ is the bare velocity of light and ^ is the bare quark mass, therefore 
a change of v t has adverse effects on masses and c(0). As the mass splittings do not change 
noticeably with the bare quark mass m or their dependence on v t has to be explained 
in some other way, which we do not yet know. 

In runs 5 and 6, the two clover term coefficients C s s w and C* w have been increased by 10% 
one at a time over their estimated values used in run 0. The observation is that an increase 
of C s s w or C* w has no noticeable effect on the value of c(0), but reduces meson masses, yet 
increases spin-spin splittings, while the effects from a 10% increase of C s s w are much larger 
than that from a similar increase in C* w . As the clover terms enter the lattice action locally 
just like mo, they are not expected to influence c(0). However, they are well-known in light 
quark calculations to make a positive contribution to meson masses []. As to the hyperfine 



1 In other words, the critical quark mass where pion becomes massless is less negative when the 
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splitting l 3 Si — l 1 ^, which comes from the spin-spin interaction between two charm quarks, 
we first note that the spatial clover terms corresponds to lattice corrections to the chro mo- 
magnetic coupling a ■ B. Hence as is expected, the spin splitting is subject to the values of 
the CiL clover coefficient. 

H. ^-tuning vs. ^-tuning 

The above observation of the influence of C s s w and Ct, w on the mass splittings is very 
important, especially since in our calculations C s s w and C* w are only (possibly very well) 
estimated. While the quenching effects will remain, some percentages of the discrepancies 
between our results and the experimental data may simply disappear when the numerical 
determination of C s s w and C* w becomes feasible in the future. While this determination (say, 
by applying the Schrodinger functional on an anisotropic lattice for heavy quarks) may be 
daunting both theoretically and computationally, this discussion of the sensitivity of our 
results to Cg W and C!, w leads us to comment on the choice of z/f-tuning over ^-tuning in this 
work. 

In the initial work |||7|, most calculations were done with the ^-tuning. If C* w and 
C* w are known numerically for the values of m and (3 that enter into our simulations, [] 
it should not matter which way is chosen to restore the relativity of the anisotropic lattice 
action, results from these two ways should agree even at finite lattice spacings up to 0(a 2 ) 
errors. If for both ^-tuning and ^-tuning, C^ w and C* w have weak or linear dependency on 
atmo as we certainly have hoped for, it would still be legitimate for the a 2 extrapolation to 
apply, so that at zero lattice spacing we would end up with the same results. Unfortunately 
in the earlier work ||[7]], ^-tuning and z/ s -tuning were not found to give the same results 
in the continuum limit, at least not by using the clover coefficients estimated from eq. (|30|) . 
Therefore at least in one of these two approaches, the dependence of C| w and C!, w on a t mo 
may be too strong to justify the a 2 extrapolation for the values of a studied here. 

Based on the classical improvement discussed in section [111 C| , we suspect that: 1. The 
two tunings would have been much more nearly consistent if the tadpole improvement was 
not based on the tree-level estimate given in eq. (|T5|)-(|16D, but instead came from the better 
classical estimate given by eq . (p~7| ) - (pT8|) . 2. If we estimate the clover coefficients by applying 
tadpole corrections to the possibly less precise estimates in eq.(|i5|)-(|l6|), which has been 
our practice so far, the results from z/ r tuning should be more reliable than the results from 
z/ s -tuning. 



clover terms are added. 

2 Or the other way around, which is less likely since in isotropic cases the numerical clover coef- 
ficients are found to be larger than the tree-level estimates, thus the splittings will be increased 
toward experimental data. 

3 Right now we ignore the £ dependence of C s s w and C* w as we have fixed the value of the 
renormalized anisotropy £. 
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Here is why. Without numerical determination of the clover coefficients, the crucial 
assumption or hope underlying the a 2 continuum extrapolation is that the clover coefficients 
depend on mass weakly or linearly. Classically for the z/ s -tuning both C* w and C| w depend 
on m c a, while for the ^-tuning only C* w depends on m c a and there is no mass dependence 
in C s s w — i/ s — 1. Furthermore, looking at the field redefinition in eq. fl4"T|) , we find it 
contradictory in estimating Cg W and C* w to use eq. (|30|) for both tunings. While this has 
been the practice so far, the clover coefficients would be effectively larger for the v s - than 
for the z/ r tuning if v t > 1 (thus the mass splittings would be larger too) and the other way 
around if v t < 1, which indeed is what is qualitatively found in ||[/J. From what we see 
in the 10% change test for the clover coefficients, the resulting discrepancy should be quite 
pronounced since v s or v t deviates from 1 by 1-12%, the range described in table |I[ 

In short, we expect the choice of tuning to be a minor issue if the clover coefficients 
are estimated in the better way described above, and it would not be an issue if C| w and 
C* w were known numerically. Most likely, on an isotropic lattice this problem of the mass 
dependence of C s s w and C^ w is only going to be more severe. In retrospect, we should still 
first tune the bare velocity of light with the clover coefficients estimated from eq.(|TTD-(|T6|). 
but once we know the tuned values of v t or v S) we should plug them into eq.(|17l)-(|T8"l) to 
get a better estimate of clover coefficients, and then use the better estimate in the following 
tunings and real runs. 



VI. CONCLUSION 

By running quenched simulations using an anisotropic action, we have been able to 
predict more reliable masses within the charmonium family than has been done in previous 
lattice calculations. The masses of both the radial n = 1 ground state and the n = 2 first 
excitation have been computed for the particles r\ c ( 1 S'o), J/ij) ( 3 Si), h c ( 1 -Pi), Xco ( 3 Po), 
and Xci ( 3 Pi)- On a finer scale, from the continuum extrapolation we have the S"-wave 
hyperfine splitting Am 1 3 Sl _ 1 i S(J of 71.8(20) MeV, the P-wave fine structure /\mi-3. Pl _izp of 
65(3) MeV, and the IP — IS" splitting Am 1 ip 1 _ 1 5 of 431(3) MeV, which agrees with other 
lattice approaches p8j , pT9[ . 



Our work shows the intrinsic benefit of an anisotropic lattice where the temporal lattice 
spacing is finer than the spatial one. At relatively low computational cost, on an anisotropic 
lattice the signals of a hadron correlator are good on more time slices. This is important 
to the calculations involving heavy quarks or excited states as they die out fast on current 
isotropic lattices. The space-time exchange symmetry, broken both on a heavy-quark ac- 
tion and (only more explicitly) on an anisotropic lattice, has been restored by tuning the 
bare parameter v t based on the dispersion relation without resorting to the "kinetic mass 
prescription" . 

While all errors given in data tables are statistical errors, the biggest errors in our results 
should be attributed to the systematic errors from quenching. Besides that, we have not 
numerically determined the two clover term coefficients C s s w and C* w , while the numerical 
tuning has been done for all other simulation inputs. It will require significant theoretical 
and computational effort to get rid of errors coming from these two sources, which is equally 
true in other lattice approaches to heavy quark systems. A feasible project in the near 
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future is to run more simulations at finer lattice spacings. By doing so we will either have 
more support for current estimations of C s s w and C* w , or we will see the breaking of the a 2 
extrapolation and thus be forced to pursue a fully numerical determination of C| w and C* w . 
Either way progress will be made. 



ACKNOWLEDGMENTS 



The numerical calculations were done on the 400 Gfiop QCDSP computer [2(J at 



Columbia University. This research was supported in part by the DOE under grant # 
DE-FG02-92ER40699. Ping Chen would like to thank Prof. Norman H. Christ for being 
Ping's Ph.D. thesis advisor. She is grateful also to all other current and former Columbia 
QCDSP members including Prof. Robert D. Mawhinney, Dr. Dong Chen, Calin Cristian, 
George Fleming, Dr. Chulwoo Jung, Dr. Adrian Kaehler, Xiaodong Liao, Guofeng Liu, 
Dr. Yubing Luo, Dr. Catalin Malureanu, Dr. Thomas Manke, Chengzhong Sui, Dr. Pav- 
los Vranas, Lingling Wu and Yuri Zhestkov. At last, this work would not have been here 
without the significant contribution from Dr. Tim Klassen. 



23 



REFERENCES 



[1] G. P. Lepage et al, Phys. Rev. D46, 4052 (1992). 

[2] C. Davies hep-ph/9710394. 

[3] H. D. Trottier, Phys. Rev. D55, 6844 (1997). 

[4] N. H. Shakespeare and H. D. Trottier hep-lat/9802038. 

[5] A. X. El-Khadra, A. S. Kronfeld, and P. B. Mackenzie, Phys. Rev. D55, 3933 (1997), 
hep-lat/9604004. 

[6] T. R. Klassen, Nucl. Phys. B Proc. Suppl. 73, 918 (1999), hep-lat/9809174. 

[7] T. R. Klassen Heavy quarks on anisotropic lattices (unpublished). 

[8] T. R. Klassen, Nucl. Phys. B533, 557 (1998), hep-lat/9803010. 

[9] M. Creutz, Phys. Rev. D21, 2308 (1980). 
[10] A. D. Kennedy and B. J. Pendleton, Phys. Lett. 156B, 393 (1985). 
[11] N. Cabibbo and E. Marinari, Phys. Lett. 119B, 387 (1982). 
[12] R. Sommer, Nucl. Phys. B411, 839 (1994), hep-lat/9310022. 

[13] R. G. Edwards, U. M. Heller, and T. R. Klassen, Nucl. Phys. B517, 377 (1998), hep- 
lat/9711003. 

[14] R. G. Edwards, U. M. Heller, and T. R. Klassen Unpublished. 

[15] M. R. Hestenes and E. Stiefel, Journal of Research of the National Bureau of Standards 

Vol. 49, No.6, December 1952 Research Paper 2379. 
[16] W. Press, B. Flannery, S. Teukolsky, and W. Vetterling Numerical Recipes in C: The 

Art of Scientific Computing. 
[17] C. Caso et al. (Particle Data Group), European Physical Journal C3, 1 (1998). 
[18] J. Simone Private communication on Fermilab results. 
[19] P. Boyle (UKQCD collaboration) hep-lat/9903017. 

[20] D. Chen et al, Nucl. Phys. Proc. Suppl. 73, 898 (1999), hep-lat/9810004. 



24 



TABLES 



p 


5.6 


5.7 


5.7 


5.9 


6.1 


L 3 • T 


8 3 • 32 


8 3 • 32 


16 3 • 64 


16 3 • 64 


16 3 • 64 


£o (£ = 2) 


1.632156 


1.654729 


1.654729 


1.690713 


1.718306 


m • a t 


0.69 


0.51 


0.51 


0.195 


0.05 


(I s 

^ sw 


2.364 


2.138 


2.138 


1.889 


1.7614 


fit 


1.429 


1.3252 


1.3252 


1.2055 


1.1431 


v s 


1 


1 


1 


1 


1 




0.92 


1.01 


1.01 


1.09 


1.12 


^ jj LandauGauge ^ 


0.7504(2) 


0.7762(2) 


0.7762(2) 


0.8091(2) 


0.8280(2) 


/ j ,-LandauGauge \ 
\ U t ) 


0.9321(1) 


0.9394(1) 


0.9394(1) 


0.9504(1) 


0.9569(1) 


c(0) 


1.012(2) 


1.000(2) 


0.991(3) 


0.984(3) 


0.984(3) 


m(15)i at (GeV) 


3.063(1) 


3.079(2) 


3.078(2) 


3.069(2) 


3.044(2) 


m(15) exp (GeV) 


3.0676(1) 


3.0676(1) 


3.0676(1) 


3.0676(1) 


3.0676(1) 


configurations 


1480 


1350 


440 


410 


588 


spatial L (fm) 


2.02 


1.66 


3.32 


2.17 


1.54 



TABLE I. This table consists of four parts. From top to bottom: 1. basic simulation inputs; 2. 
mean links measured in Landau gauge used to estimate the clover coefficients; 3. tuning accuracy of 
the effective velocity of light c(0) and the IS spin average meson mass; 4. number of measurements 
and the physical spatial extension. All the errors on masses, m(lS) in this table and more masses 
in the following tables, are purely statistical errors. 



p 


5.6 


5.7 


5.9 


6.1 


z 


2 


2 


2 


2 


Co 


1.632156 


1.654729 


1.690713 


1.718306 




1.982(10) 


2.413(6) 


3.690(11) 


5.207(29) 


L (fm) if ± = 8 


2.018 


1.658 


1.084 


0.768 


L (fm) if ^ = 16 

v ' a s 


4.036 


3.315 


2.168 


1.536 


1/at in GeV 


1.564 


1.905 


2.913 


4.110 



At a given f3 and true anisotropy £, the corresponding bare anisotropy £0 and the 

Both were determined to 1% accuracy. On the 
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TABLE II 

Sommer scale tqn/o-s are quoted from 
last line 1/at is given in GeV, numbers which are used in this work to express the charmonium 
mass spectrum in physical units. 
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r 


25+1 L 


J PC uu 


cc 


75 


'So 


_+ 7T 


Vc 


7s 


3 S, 


1~ P 


j/y 


Is Is, 


'Pi 


1+- bl 


h c 


1 


3 Po 


0++ a 


XcO 


757s 


3 Pi 


1++ ai 


Xcl 


TABLE III. Meson states created by local and relativistic operators of the form ipTifj, labelled 


in spectroscopic notation and by the particles they are identified with in the system of 1: 


ight flavors 


uu and the charmonium family cc. 






P 


5.6 


5.7 5.7 5.9 


6.1 


L 3 • T 


8 3 • 32 


8 3 • 32 16 3 • 64 16 3 • 64 


16 3 • 64 


Therm sweeps 


2500 


2500 2500 5000 


5000 


Landau gauge fixing 


all elements of the traceless antihermitian part < 10' 


-5 


updating 


100 


100 100 500 


500 


Landau link config 


80 


15 15 10 


10 


Spectrum updating 


100 


100 100 200 


400 


sink & source 




local sink, three box sources of different sizes 




source (N x ■ N y ■ N z ) s 


2-2-2 


2-2-2 2-2-2 3-3-4 


4-4-4 


source (N x ■ N y ■ N z ) m 


3-3-3 


4-4-4 4-4-4 5-6-6 


7-8-8 


source (N x ■ N y ■ N z )i 


6-6-6 


6-6-6 6-6-6 11-11-11 


13-13-13 


Coulomb gauge fixing 




sum of square of all 18 matrix elements < 10 -9 




CG stopping end 


1.0e-8 


1.0e-8 1.0e-8 1.0e-8 


1.0e-8 


CG iterations 


42(4) 


43(3) 49(4) 64(4) 


85(5) 


Mean plaq (D) s 


0.41318(7) 


0.43925(6) 0.43918(3) 0.47723(2) 


0.50295(2) 


Mean plaq (D) t 


0.65693(3) 


0.67401(3) 0.67398(1) 0.69923(1) 


0.71722(1) 


TABLE IV. Simulation paramet 


ers detailed in six parts from top to bottom: 1. 


inputs to 


identify each simulation; 2. heatbath ; 


sweeps for thermalization; 3. Landau gauge fixing 


; condition, 



heatbath updating sweeps, and number of measurements of mean links; 4. heatbath updating 
sweeps, sink and source types, small-size box source, medium-size box source, large-size box source, 
and Coulomb gauge fixing condition for the mass spectrum; 5. stopping condition and number 
of iterations of the conjugate gradient (CG) in the calculation of quark propagators; 6. measured 
average values of spatial and temporal plaquettes. 
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5.6 


5.7 




5.9 


6.1 


Nature 


L 3 ■ T 


8 3 -32 


8 3 • 32 


16 3 • 64 


16 3 • 64 


16 3 • 64 


N/A 


m i^s (GeV) 


3.001(1) 


3.023(2) 


3.021(2) 


3.012(1) 


2.992(2) 


2.9798(2) 


m 2hS ( GeV ) 


3.61(3) 


3.70(6) 


3.65(3) 


3.70(2) 


3.66(3) 


3.594(5) 


m 1 3 Sl (GeV) 


3.084(1) 


3.099(2) 


3.098(2) 


3.090(1) 


3.062(2) 


3.09688(4) 


m 2 3 Sl (GeV) 


3.65(3) 


3.64(3) 


3.78(2) 


3.75(4) 


3.73(3) 


3.68600(9) 


m 1 ip 1 (GeV) 


3.523(3) 


3.526(5) 


3.518(5) 


3.517(5) 


3.498(12) 


3.5261(2) 


m 2 i Pl (GeV) 


4.12(5) 


3.92(10) 


4.02(5) 


4.09(5) 


3.86(10) 


N/A 


m l3Po (GeV) 


3.499(2) 


3.481(21) 


3.480(2) 


3.465(3) 


3.421(4) 


3.417(3) 


m 2 3 Po (GeV) 


3.82(5) 


3.79(8) 


4.05(5) 


4.03(4) 


4.02(5) 


N/A 


m 1 3p 1 (GeV) 


3.530(2) 


3.523(6) 


3.518(4) 


3.519(2) 


3.472(6) 


3.5105(1) 


m 2 3 Pl (GeV) 


3.92(5) 


3.85(11) 


4.08(6) 


3.99(5) 


4.00(8) 


N/A 



TABLE V. Charmonium spectrum measured at four values of f3, compared with their observed 
values in nature. As explained in the text, the errors from scale setting are not included, which 
also applies to the tables that follow. 









5.6 


5.7 




5.9 


6.1 


a -» 


Nature 


L 3 -T 






8 3 • 32 


8 3 • 32 


16 3 • 64 


16 3 • 64 


16 3 • 64 


N/A 


N/A 


Am l3Sl . 




(MeV) 


84.7(4) 


75.9(5) 


77.2(6) 


76.4(6) 


73.3(6) 


71.8(20) 


117.1(2) 


Am^ p 1 . 


-i 3 p 


(MeV) 


33(1) 


37(3) 


43(3) 


59(1) 


58(2) 


65(3) 


93(3) 


Am 1 ip 1 . 


-is (MeV) 


473(7) 


449(6) 


454(6) 


442(5) 


439(5) 


431(3) 


458.5(2) 






(MeV) 


611(28) 


675(57) 


625(27) 


691(20) 


666(34) 


693(19) 


614(5) 


Am 2 3 Sl . 


-135! 


(MeV) 


563(32) 


546(29) 


677(20) 


662(35) 


665(34) 


696(40) 


589.1(1) 


Am 2 i Pl _ 




(MeV) 


591(48) 


391(103) 


503(48) 


575(52) 


358(102) 


412(92) 


N/A 


Am 23 p „ 


-l 3 Po 


(MeV) 


321(52) 


307(76) 


567(50) 


566(41) 


603(53) 


667(75) 


N/A 


Am 2 3 Pl . 


-13 Pi 


(MeV) 


385(52) 


327(105) 


559(55) 


476(53) 


532(78) 


549(71) 


N/A 



TABLE VI. Mass splittings in the charmonium spectrum measured at four values of (3, their 
extrapolated values to the continuum limit, and the observed values in Nature. 
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r ltteo. 


— — — 

fitting 


2 

X per 






goodness 




dropped 


J „ -C 

a. o.i. 


bin 


quantities 


range 


Q.O.I. 


value 


oi nt 


eigenvalues 


leit 


size 


c(p = 0) 


6-13 


2.8(5) 


0.003 


0.029 


9 


9 


2 


mis, Am l35l _ii5 , Amii Pl _is 


6-13 


1.6(5) 


U.ll 


1.1 


Q 

8 


1 n 
10 


2 


/\ 777,1 '(n i q o 


7-13 


7(5) 


0.73 


7.3 


o 


10 


2 


m^so, Am 2 is _iis 


4-13 


0.8(4) 


0.68 


10.1 


8 


15 


2 


mi3 5l , Am 2 3 Sl . l3Sl 


5-13 


1.3(4) 


0.22 


3.1 


6 


14 


2 


mii Pl , Am, 2 i Pl _ii Pl 


2-13 


1.5(3) 


0.06 


1.2 


8 


21 


2 


mi3 Po , Am, 2 3 Po _i3 Po 


3-13 


1.0(3) 


0.42 


9.5 


3 


23 


2 


mi3 Pl , Am 2 3 Pl _i3 Pl 


3-13 


1 2f3l 


0.26 


6.0 


3 


23 


2 


TABLE VII. Fitting details at (3 = 


(J.U UII dll < 


S 3 • 32 lattice. See the text for details on 


UIlc 


fitting ansatz, procedure and criteria. 














Fitted 


fitting 


2 

X per 


n 
H 


goodness 


dropped 


H n f 


bin 


quantities 


range 


U.O.I. 


va 1 n p 


of fit 


eigenvalues 


left 


size 


c(d = 0) 


8-15 


U.y(4 ) 


n 58 


7.5 


5 


1 3 


Z 


mis, Am^Si-iiSo; Amii Pl _is 


7-15 


1.3(4) 


99 


3.1 


7 


14 


2 


Am l3Pl _ l3Po 


7-15 


1 2(4) 


94 


3.4 





14 


2 


m i^s i Am 2 is _iis 


8-15 


1.5(4) 


u. 10 


1.6 


5 


1 9 


2 


mi3 Sl , Am 2 3 5l _ l35l 


8-15 


1.2(4) 


0.26 


3.8 


2 


15 


2 


mii Pl , Am 2 i Pl _ii Pl 


5-15 


1.0(3) 


0.44 


8.8 


6 


20 


2 


m l3Po , Am 2 3 Po _i3 Po 


5-15 


1.6(3) 


0.04 


0.9 


5 


21 


2 


m l3Pl , Am 2 3 Pl _i3 Pl 


5-15 


0.9(3) 


0.63 


13 


6 


20 


2 


TABLE VIII. 


Fitting 


CLcLdllo djL p 


' = 5.7 


on an 8 3 • 32 lattice. 






Fitted 


fitting 


2 

X per 


Q 


goodness 


dropped 


d.o.f. 


bin 


quantities 


range 


U.O.I. 


value 


of fit 


eigenvalues 


left 


size 


c(p = 0) 


6-16 


2.5(3 ) 


0.001 


0.02 


13 


14 


2 


mis, Ampsi-iiso' Amii Pl _is 


5-15 


1.4(3) 


0.14 


2.5 


9 


18 


2 


Ami3 Pl _i3 Po 


5-15 


0.7(4) 


0.85 


14 


2 


16 


2 


"iiiSoi Am 2 is -iis 


4-15 


1.0(3) 


0.42 


8.0 


10 


19 


2 


mi3 5l , Am 2 3 Sl _i3 5l 


3-15 


1.1(3) 


0.38 


7.2 


13 


19 


2 


mii Pl , Am 2 i Pl _ii Pl 


3-15 


1.2(3) 


0.19 


4.7 


8 


24 


2 


mi3 Po , Am 2 3 P() _i3 P() 


3-15 


1.2(2) 


0.19 


5.9 





32 


2 


mi3 Pl , Am 2 3 Pl _i3 Pl 


3-15 


1.1(3) 


0.31 


7.5 


8 


24 


2 



TABLE IX. Fitting details at (3 = 5.7 on a 16 3 • 64 lattice. 
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r ltted 


— — 

fitting 


2 

X per 






goodness 




dropped 


^1 ^ -p 

Q.O.I. 


bin 


quantities 


range 


Q.O.I. 


value 


01 nt 


eigenvalues 


leit 


size 


c(p = 0) 


9-24 


1.8(3) 


0.009 


0.21 


19 


23 


2 


mis, Am l35l _ii5 , Am^p^s 


o-z4 


1.2(3) 


0.1 f 


4.y 


16 


on 
29 


z 


Ami3p, _i3r. 


5-24 


1.3(2) 


0.12 


4.4 





36 


2 


m i'So. Amaiso-iiso 


6-24 


1.2(2) 


0.18 


5.8 


17 


33 


2 


m l3Sl , Am 2 3 Sl . l3Sl 


7-24 


1.5(2) 


0.04 


1.2 


15 


32 


2 


m^pj, Am 2 i Pl _ii Pl 


5-24 


1.2(2) 


0.22 


8.2 


15 


38 


2 


m l3Po , Am 2 3 Po _i3 P() 


5-24 


1.4(2) 


0.05 


2.3 


11 


42 


2 


m l3Pl , Am 2 3 Pl _i3 Pl 


6-24 


1.0(2) 


0.54 


23 


7 


43 


2 



TABLE X. Fitting details at (3 = 5.9 on a 16 3 • 64 lattice. 



Fitted 


fitting 


X 2 Per 


Q 


goodness 


dropped 


d.o.f 


quantities 


range 


d.o.f. 


value 


of fit 


eigenvalues 


left 


c(p = 0) 


17-32 


1.3(3) 


0.15 


4.1 


14 


28 


mis, Am^^i^, Am^p^s 


13-32 


1.2(2) 


0.16 


5.8 


17 


37 


Am l3Pl _ l3Po 


12-32 


1.2(3) 


0.18 


5.5 


7 


31 


n*liS > Am 2i5 -l 1 5 


14-32 


1.0(2) 


0.53 


20 


13 


37 


m l35l , Am 2 3 Sl _ 1 3 5l 


14-32 


0.9(2) 


0.60 


22 


13 


37 


m!i Pl , Am 2 i Pl _ii Pl 


11-32 


2.1(2) 


10- 5 


4 x 10~ 3 


9 


50 


m l3Po , Am 2 3 Po _ 1 3 Po 


10-32 


1.5(2) 


0.02 


0.7 


16 


46 


m l3Pl , Am 23Pl _i3 Pl 


10-32 


1.3(2) 


0.08 


3.4 


17 


45 



TABLE XI. Fitting details at (3 = 6.1 on a 16 3 • 64 lattice. 



Meson Sz momenta 


c(0) 


fitting 


X 2 Per 


Q 


goodness 


dropped 


d.o.f. 


bin 






range 


d.o.f. 


value 


of fit 


eigenvalues 


left 


size 


1%, (1,0,0) (2,0,0) 


0.984(3) 


17-32 


1.3(3) 


0.15 


4.1 


14 


28 


2 


1%, (1,1,0) (2,2,0) 


0.979(14) 


14-32 


1.5(3) 


0.06 


1.6 


24 


27 


2 


l 3 ^, (1,0,0) (2,0,0) 


0.986(3) 


17-32 


1.3(2) 


0.08 


3.2 


1 


41 


2 


1 3 P , (1,0,0) (2,0,0) 


0.973(17) 


17-32 


1.2(2) 


0.16 


6.9 





42 


2 



TABLE XII. Dispersion relation checked at (3 = 6.1 on a 16 3 • 64 lattice. As an example on 
how to read this table, the first row gives the effective velocity of light c(p = 0), fitted from the 
pseudo-scalar meson l 1 ^ with momenta 0, ^(1, 0, 0) and ^-(2, 0, 0). Not listed here due to failed 
fitting are: meson l 3 Pi and l 1 -/^ of the (1,0,0) momentum series, and all mesons but l 1 5o of the 
(1, 1,0) momentum series. 
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ID 




Inputs 






Outputs 


(masses in 


MeV) 




#of 




m 


v-t 






c(0) 


15 


AS 


AF 


ASP 


conf 





0.51 


1.01 


2.178 


1.3396 


0.996(2) 


3070(2) 


79.1(6) 


43(3) 


455(8) 


980 


1 


0.5355 


1.01 


2.178 


1.3396 


1.000(6) 


3125(3) 


78.0(9) 


40(5) 


467(13) 


190 


2 


0.4845 


1.01 


2.178 


1.3396 


1.001(4) 


3017(3) 


82.9(9) 


40(6) 


463(13) 


200 


3 


0.51 


1.06 


2.178 


1.3396 


0.974(5) 


2991(2) 


75.9(9) 


30(5) 


470(30) 


200 


4 


0.51 


0.96 


2.178 


1.3396 


1.035(5) 


3156(3) 


82.4(10) 


39(5) 


440(15) 


230 


5 


0.51 


1.01 


2.396 


1.3396 


1.001(3) 


3033(2) 


92.8(7) 


42(4) 


463(8) 


730 


6 


0.51 


1.01 


2.178 


1.4736 


1.000(2) 


3042(2) 


83.3(6) 


45(5) 


462(8) 


470 



TABLE XIII. At [3 = 5.7, £ = 2 on an 8 3 ■ 32 lattice, we check the effects of the inputs on the 
effective velocity of light c(0), the IS spin average meson mass m(lS), the S-wave mass splitting 
Am 1 35 1 _ 1 i5 , the P-wave mass splitting Am 1 3p 1 _ 1 3p , and the S — P splitting Am 1 ip 1 _ 1 5- The 
right-most column gives the number of configurations. Inputs that are not shown are exactly the 
same as those listed for the (3 = 5.7 runs in table | and table IV. 
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FIG. 1. Measured at four values of lattice spacings, the charmonium masses are plotted on 
the top of experimental values (shown as horizontal lines wherever available). 



31 



1 3 S,-1 1 S 



1 3 P,-1 3 P 




1 P, - 1S 



2 1 S n -1 1 S„ 





0.00 0.01 0.02 0.03 , 0.04 0.05 0.06 

a s 2 (fm 2 ) 



0.00 0.01 0.02 0.03 , 0.04 0.05 0.06 

a s 2 (fm 2 ) 



FIG. 2. The continuum a s extrapolation of mass differences. The four squares correspond to 
the four f3 values, 5.6, 5.7, 5.9 and 6.1, at fixed renormalized anisotropy £ = 2. The vertical dashed 
line is there to emphasize the continuum limit. Wherever available, a circle on the dashed line is 
the observed mass difference in Nature. 
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FIG. 4. The — IS" splitting from several approaches is presented in the form of the ratio 



Am, 



4 58 5MeV S • ^ ee sec ti° n VP for more details. The vertical dashed line is to emphasize 
the continuum limit. The solid line denotes the a? continuum extrapolation of our work. 
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FIG. 5. The l 3 S\ — l 1 5o splitting from several approaches. The vertical dashed line is to 
emphasize the continuum limit. The solid line denotes the a 2 continuum extrapolation of our 
work. 
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FIG. 6. The l 3 Pi — l 3 Po splitting from several approaches. The vertical dashed line is to 
emphasize the continuum limit. The solid line denotes the a 2 continuum extrapolation of our 
work. 
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